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Abstract. 

An efficient numerical quadrature is proposed for the approximate calculation of the 
potential energy in the context of pseudo potential electronic structure calculations with 
Daubechies wavelet and scaling function basis sets. Our quadrature is also applicable in 
the case of adaptive spatial resolution. Our theoretical error estimates are confirmed by 
numerical test calculations of the ground state energy and wave function of the harmonic 
oscillator in one dimension with and without adaptive resolution. As a byproduct we derive 
a filter, which, upon application on the scaling function coefficients of a smooth function, 
renders the approximate grid values of this function. This also allows for a fast calculation 
of the charge density from the wave function. 

1. Introduction. 

Gaussians and plane waves are at present the most popular basis sets for density 
functional electronic structure calculations. Wavelets are a promising new basis set 
that combines most of the theoretical advantages of these two basis sets. They can 
form a systematic orthogonal basis set that allows for adaptivity, the basis functions 
being localized both in real (compact support) and in Fourier space. 

The first attempts to use wavelets in the electronic structure calculations appeared 
more than 10 years ago. The first papers we are aware of used the Mexican hat 
wavelet [1], [2] and the Meyer wavelet [3]. However, these wavelet families were soon 
abandoned because they do not have compact support. Daubechies [4] wavelets were 
then investigated in a series of publications [5]- [8]. This basis is orthogonal and has 
the property of having the highest number of vanishing moments for the given support 
width, thus combining locality and approximating power. It is also localized in the 



2 



momentum space. To use the Daubechies wavelets in the variational Galerkin method 
for the Schroedinger equation, one has to compute the matrix elements of the kinetic 
and potential energy operators. The algorithm for the kinetic part is straightforward 
[9]. 

The main difficulty is the calculation of the potential energy matrix elements [10]- 
[12]. They were computed by expanding the potential in terms of scaling functions, 
too, and then a convolution was performed with the matrix of products of three scaling 
functions. We will call this "the triple product method": for details, see the beginning 
of Section 4. It requires a lot of computer resources; this motivated alternative 
approaches. The collocation approach [13] is not well suited for the Daubechies' 
wavelets since it spoils the favorable convergence rate of variational schemes. Another 
approach [14] involved designing a quadrature for the product of two scaling functions 
and a smooth function. It decreased the amount of computations in comparison to 
the triple-product method, but not sufficiently. 

It follows from the above considerations that the Daubechies basis set can only 
be useful for electronic structure calculations if one has a better algorithm for the 
calculation of the potential energy than those available at the moment. Such an 
algorithm will be proposed in the present paper. 

Due to the above listed problems with the Daubechies family, interest focused 
recently on the interpolating Deslarier-Dubuc [15] family [16]-[19]. Because the scaling 
functions of this family are interpolating (cardinal), the collocation approximation is 
much more accurate for them than for the Daubechies family. An even more accurate 
approximation for the potential energy is based on a relation with the analytically 
known overlap matrix elements [19]. 

The major disadvantage of the Deslarier-Dubuc wavelets is that they are not 
orthogonal. For very large systems, the dominating term in independent particle 
electronic structure calculations is the orthogonalization of the one particle orbitals. 
The prefactor for this dominating cubic term is much smaller if an orthogonal basis set 
is used compared to a non-orthogonal basis set. Orthogonal wavelets are in addition 
also interesting candidates for the implementation of linear scaling algorithms [20]. 

Alpert [21], [22] polynomial multiwavelets overcome the above mentioned disad- 
vantages. The potential energy can be calculated easily and they are orthogonal. It 
seems that they are the ideal basis set for all electron electronic structure calcula- 
tions and impressive results have been reported [23]- [25]. The Chui-Lian [26] family 
has also been used in the same context [27]. Since multiwavelets can represent dis- 
continuous functions, they are well suited to represent the electron-nucleus cusp in 
all electron calculations. However, if one uses the Bachelet-Hamann-Schulter [28] or 
Gaussian [29] pseudopotentials, the wavefunctions and the potential are smooth and 
this property is not needed. Then one may prefer the Daubechies wavelets due to their 
simplicity and small support length. For this reason we explore in this paper the use 
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of Daubechies wavelets for pseudopotential electronic structure calculations, together 
with our novel quadrature scheme. A topic that is closely related to the problem of 
integrating the potential energy is the problem of finding quadrature schemes for the 
product of a scaling function(wavelet) and a smooth function [14], [30]- [33]. 

The quadrature scheme to be presented in this paper is aimed at applications in in- 
dependent particle schemes such as density functional theory where one 3-dimensional 
single particle orbital is provided for each electron in the system. For completeness, 
we will mention that wavelets have also been explored as a basis set for high dimen- 
sional many-electron wavefunctions [34]- [36]. Another approach is to use the so-called 
Weyl-Heisenberg wavelets that do not have multiresolution properties but are related 
to the structure of the phase space [37]- [40]. Although these two approaches are 
promising, any treatment of correlation entails an important increase of the numer- 
ical effort and such approaches will not allow to treat systems with several hundred 
atoms in the near future. 

In the present paper we apply one-scaling-function quadratures to the numerical 
calculation of the potential energy matrix element (1) between two smooth functions. 
For this purpose we developed an algorithm for the reconstruction of grid values of 
a smooth function from its scaling function expansion coefficients. This technique 
might also be used in other contexts such as in speech reconstruction. Our scheme 
also provides a way to calculate the density from the wave function expressed in 
scaling functions. The extension of our potential energy quadrature onto the case of 
adaptive spatial resolution is then described. 

The paper is organized as follows: 

In Section 2 we briefly recall the definition and properties of the Daubechies 
wavelet family. 

In Section 3 we recall the quadrature of [14], [30] for the product of a scaling 
function and a smooth function. 

Using that, in Section 4 we construct the quadrature for the potential energy 
functional, i.e., for the product of the form 

< smooth function \ potential \ another smooth function > . (1) 

In subsection 4.1 we derive the quadrature and estimate its error, which behaves 
essentially as a square of the error for the wavelet expansion of the smooth functions 
involved. Then in subsections 4.2-4.3 we prove that the quadrature is exact if the 
potential and one of the functions in (1) are polynomials and another function is a 
scaling function. This confirms the previous estimation of error and extends it to the 
case when only one of the functions in (1) is smooth. 

In Section 5 we extend our method onto the case of adaptive resolution, repeating 
the procedure of the previous Section. The errors are again estimated. It is shown 
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that now the main source of error is the boundary between regions with different 
resolution. 

In Section 6 we modify our adaptive quadrature so that in the regions with con- 
stant resolution it reduces to the non-adaptive one. This reduces the computational 
burden. The price is that now we need to minimize the Rayleigh-Ritz (RR) func- 
tional in the space of smooth functions, i.e., smoothen the cusp of its gradient at the 
boundary between regions with different resolution. 

In Section 7 we present a way to compute the density corresponding to a wave 
function expressed in scaling functions. It uses the approximate wave function values 
derived by the method described in Appendix A. This way is fast, but it reproduces 
multipole moments of the density only approximately (although with good precision). 

Finally, in Section 8 we apply our methods to the calculation of the ground state 
energy and wave function of the harmonic oscillator, both with and without adap- 
tivity. The results are then compared with those obtained with the RR functional in 
which the potential energy was calculated exactly. We consider both the least asym- 
metric and extremal phase Daubechies wavelets and argue that the least asymmetric 
family is preferable. 

2. The orthogonal wavelets. 

In this work we use the Daubechies [4] scaling functions <f>(x) and wavelets ij)(x), in 
the dilated and shifted form: 

<j) k (x) = 2 k/2 <j)(2 k x - i); $(x) = 2 k/2 ^(2 k x - i). 

where i, k are integer. Sometimes we will also employ the intermediate notation: 

<j>i(x) = (f>(x - i); (j) k (x) = 2 k/2 (j)(2 k x). 

Our conclusions are the same for the least asymmetric and extremal phase [4] Daubechies 
wavelets, but the least asymmetric ones behave better in the examples we considered. 
Thus we will use the least asymmetric family for the illustration. The graph of the 
scaling function of this family of the order 8 is given on Fig. I. 

Since the Daubechies scaling functions and vectors are orthogonal, it is natural to 
use the bra and ket notation for them. Then, their orthogonality conditions can be 
written as 

< $1$ >= Sij] <0 J fc |^f>=O, k' > k; < tftyf >=6 ij 6 kk >. 
One can define the following pair of sequences of spaces: 



V k = span{ | <p k > } ; W k = span{ \ i\) k > } 



(2) 
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where Vk n V\4 = {0}. It turns out that for any k > the space of square integrable 
functions C 2 can be decomposed into the following infinite direct sum: 

C 2 = V k © W k © W fc+ i © . . . . (3) 

In this paper we will need the refinement relations: 

l^t 1 >= E h M i+j >; l^" 1 >= E^lfe > • (4) 

They are just a unitary transformation to a new basis; for the spaces (2) it means 
that 

The inverse of (4) is called the forward wavelet transformation, but we do not need 
its explicit form in the present text. 

The decomposition (3) can be reformulated in the following way: any square 
integrable function f(x) can be uniquely expanded as 

oo 

/(*) = E + */(*); = E E <WK X )- ( 6 ) 

j q=k j 

The Daubechies-2m scaling functions are a basis of degree m — 1. This means 
that every polynomial of degree less than m is contained in V^. Therefore [41], the 
tail part of the series (6) behaves as 

5f( x ) = 0{h m ); h = T k (7) 

with respect to the C 2 norm. The O notation corresponds to the limit of k — > oo. 
For the polynomials of degree less than m, Sf(x) = 0. 

The methods to be presented can be extended to other wavelet families, not nec- 
essarily orthogonal, but in our opinion, the Daubechies family is the optimal choice 
for the electronic structure calculations. 

3. Quadrature for orthogonal wavelets. 

A wavelet quadrature is determined by the set of coefficients wi such that for a smooth 
function G(x), the integral / G(x)(f)^(x)dx is approximated by VhJ2iWiG , f + , n where 
G) = G{2- k x). The y/h factor comes from the normalization of the scaling function. 
The degree of accuracy of a quadrature formula is M if it yields the exact result 
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for every polynomial of degree less than or equal to M. This is equivalent to the 
condition [14] 

5>,Z a = M s ; M s = I 'y s cf>(y)dy (8) 
i J 

where M s are the scaling function moments. 

If the quadrature filter is of degree M and G(x) belongs to C M+1 then [14], [30] 

[ G(x)<f> k r (x)dx = VhJ2wiGl r + 0(h M+3 / 2 ) = VhJ2w s -rG k s + 

J l s 

+ 0{h M+ *' 2 )- h = 2- k . (9) 

In this paper we will usually work with the uniform quadrature (of degree M) set 
forth in [14]: 

M 

Wl = Y,PlrM r . (10) 

r=0 

It will be supposed that M = 2m — 1 for the Daubechies-2m wavelets, and I = 
1 — m..m. In (10), the Lagrange polynomials of degree M are used: 



l-j f^L r\ dy r 



j=l—m J r=0 



The nonzero values of the filters wi for m = 3.. 6 are shown at the Table 1. 

Fig. 1 also contains the values of w t for m = 4 and M — 7, compared with the 
graph of the corresponding scaling function (least asymmetric Daubechies-8). One 
sees that the filter values are close to the scaling function values at integer points. 



4. The potential energy functional for orthogonal 
wavelets. 

An important step in the electronic structure calculation is to find the eigenspectrum 
of a one-particle Schroedinger equation. Although our approach can be extended to 
the three-dimensional systems, in this paper we will consider only one-dimensional 
ones: 

HV(x) = EV(x); H = --S-2+V(x) =f + U; [dxV 2 (x) = l. (11) 
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Table 1. The values of the filters wi for the Daubechies-2m least asymmetric filters; 
m = 3..6. 



1 



Daubcchies-6 



Daubechies-8 



Daubechics-10 



Daubcchics-12 



-5 
-4 
-3 

-2 0.0858797754503928 

-1 1.0472376804223309 

-0.1886782932535312 

1 0.0795781221430145 

2 -0.0288721312776034 

3 0.0048548465153963 
4 

5 
6 



0.0026299127476935 
-0.0377927339236569 
0.0755988357512099 
0.9999560903030736 
-0.0794124676160406 
0.0451427040622791 
-0.0069875964135745 
0.0008652550890159 



0.0003712028220936 
-0.0046529756260417 
0.0306436002784248 
-0.1207447752890374 
0.1338108260452157 
0.9123169219278740 
0.0109419516584456 
0.0393078583967683 
-0.0022599250999316 
0.0002653148861886 



0.0000754232174770 
-0.0011760498174610 
0.0104347966396891 
-0.0340901829704789 
-0.0067678682684262 
1.0005931732054807 
0.0041859363010669 
0.0351468153360141 
-0.0096794739531791 
0.0015648660417616 
-0.0003139771845937 
0.0000265414526497 



This equation can be solved by representing a trial wavefunction as a linear combina- 
tion of finite elements [42] . In the case considered in this paper, these finite elements 
are the Daubechies scaling functions (wavelets). Thus we select the trial function in 
the form 

<M*) = E ( 12 ) 

i 

which can also be represented by the ket vector \c k >= X^L-jv c i|<^ > with N = 
L/h = 2~ k L. In this paper we use nonperiodic boundary conditions, although we 
could use the periodic ones too. Then if the Rayleigh-Ritz functional is defined as 

n < c k \H\c k > T({c k }) + U({c k }) _ 
R (M) = <c u\ ck> = T^F^ > ( 13 ) 

T{{c k }) = < c k \f\c k >; U{{c k }) =< c k \U\c k > (14) 

the variational ground state energy and wave function expansion coefficients have the 
form 

E = mmR({c k }); \c k >= &rgmmR{{c k }). (15) 

The other wavefunctions are obtained in the same way using the Lagrange multipliers. 

Since the Daubechies-2m scaling functions have degree m — 1 (see Section 2), the 
ground state energy error behaves [42] as 

6E = E - E ex = 0(h 2m ~ 2 ) (16) 
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where E ex is the exact ground state energy of (11). 
The kinetic energy has the form [43] 

where ai is the kinetic energy filter. 

The potential energy (14) can be written in the conventional form as 



U = J ^/(rr)!/ (a;) ^/ (re) <ir. 



For the calculation of the excited states we have to evaluate the more general potential 
energy bilinear form: 

U = J $ I (x)V(x)V I (x)dx (17) 
where the wave function <&j(x) has the same form (12): 

*,(*) = (18) 

i 

In the electronic structure calculations the exact analytical form of the potential 
V{x) is not known. Any approximation of the potential should invoke energy errors 
that alter the Rayleigh-Ritz functional (13) as little as possible; otherwise the min- 
imization with an approximate functional will not converge to the true minimum. 
We choose the following quantitative criterion for this: the error in the potential 
that arises from approximation should be small compared to (16). When the grid 
parameter h is small, there is a convenient way to ensure that: we require that the 
approximation error behaves as h 2m . Note that, e.g., the collocation approximation 
[13] does not fulfill the latter requirement. 

The most natural way to approximate the potential is to expand it in the inter- 
polating scaling functions ^(x) [15], [43] of degree 2m: 

V{x) = , EV^(x) + SV(x)- 1 5V(x) = 0(h 2m ); (19) 

i 

0f(x) = tf(z/h-i); V^ = V(ih); h = 2~ k . 

where the error estimate (19) holds for sufficiently smooth potentials. Substituting 
the expansions (12), (18), (19) into (17), we get: 

U = Y,<%Vj*i j tiWYW^ (2°) 

ijl ' ijl 

+ 0(h 2m ); I rs = [ (f) r {x)(t) I {x)(t) s {x)dx (21) 
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where (21) is the triple product matrix mentioned in the Introduction. Unfortunately, 
the formula (20) requires too many floating point operations (flops) per grid point. 
In the simplest case when sf coincide with cf, the leading term in that number is 
3/2N 2 in one dimension and 2N 4 in three dimensions. N is the dimension of the 
matrix (21), and eigenvalue decomposition of the matrix (21) is used in the 3D case. 
For the Daubehies-2m scaling functions, N = 2m — 2, so the leading term in m of 
the number of flops per gridpoint is 6m 2 in 1 dimension and 32m 4 in 3 dimensions. 
This is unacceptable for the realistic values of m of the order of 10. 

Therefore we need some approximation that would decrease the amount of com- 
putations, but still invoke the error that asymptotically behaves as 0(h 2m ), at most. 
Such an approximation will be described in the next subsection. 



4-1 The quadrature for the case of smooth wave functions. 

Let us consider two smooth functions 3>(x), *&(x) together with their approximate 
wavelet expansion $/(x), ^i(x) (12), (18). The expansion coefficients c\ and are 
obviously given by 



c\ = J ti(x)*(x)dx; sj = J <l>*(x)*(x)dx. (22) 
Then, according to (6), (7), 

$(ar) = $j(ar) +5$(x); = V z (x) + 8V(x); (23) 

SHx) = EEOfW = W; 5*(x)='£'Eb?tf(x)=0(h m ).(U) 

k'>k i k'>k i 

From (23) we get 

$ z (x) = - 6$(x); ^i(x) = V(x) - 6V(x). 
Plugging this expression into (17), we can get after some easy calculations: 

U = U A + 5U! + 5U 2 ; Ua = J <Z>i(x)V(x)V(x)dx; (25) 

5Ui = J 5®(x)V(x)5V(x)dx; 5U 2 = - J <S>{x)V{x)5t>{x)dx. (26) 

In the next Section we will prove that 

8U X + 5U 2 = 0{h 2m ) (27) 

from the local properties of scaling functions. In this Section we will show for the two 
error terms separately that 

5U 1 = 0(h 2m ); 5U 2 = 0{h 2m ). (28) 
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The first estimate of (28) follows from (24) and the fact that individual wavelet 
coefficients d\' , b k> from the tail parts behave as 2~ fe,(m+1 / 2 ) [30]. If the potential V(x) 
is bounded then it follows from the Schwarz inequality: 

J 5&(x)V(x)5V(x)dx < max\V\ J \5$(x)\\5V(x)\dx < max\V\\\5$\\\\5^\\. 

To prove the second part of (28), let us define the function H(x) = $(x)V(x). It 
can also be interpolated by scaling functions: 

H(x) = Hj(x)+5H(x); H T (x) = £ h^x)- 

i 

h k = U k {x)H{x)dx ] 6H(x)=5252tftf(x)=0(h m ). 

J k'>k i 

Then, / H I (x)8^{x)dx = 0, and 

5U 2 = - J 5H{x)5^(x)dx = 0{h 2m ). 
We are left with the approximation (25). Substituting (18) into it, we get 

Ua = Y. c 1 I <t> k i{x)G{x)dx (29) 

i J 

where G(x) = V(x)^>(x). The scalar products in (29) can now be determined by 
applying the scaling function quadrature (9): 

UaVS/t] = J2 c i\^J2 w s-iG k s + 0(h 2m+1 / 2 )] =Y J PsG k s +0{h 2m ) (30) 

i L s - 1 s 

where 

G k q = G(hq) = V(hq)V(hq) = V k V k 
are the grid values of the function G(x). As seen from (30), the quantity 



w. 



plays the role of the quadrature for the whole function $/. 

Now we need to find the grid values of the smooth function G(x). We suppose that 
we know the grid values V k of the smooth potential. Then it remains to reconstruct 
the grid values ty k of the unknown wave function from the known coefficients s k . 
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The easiest way to obtain the grid values would be by using (12): 

i 

Since the Daubechies scaling functions are not very regular, the value of a scaling 
function at a real space grid point does not very well represent the behavior of the 
scaling function in a small interval around this grid point. As a consequence of this 
and in accordance with (7), the error in (31) behaves like 0(h m ). which is much 
worse than (27) and we therefore discard this possibility. We will instead introduce 
some smoothed grid values which better represent the average behavior of 
In appendix A it will be shown that there exists a filter Wj such that the smoothed 
grid values can be obtained by a convolution from the scaling function expansion 
coefficients 

*: = 7s?«^-7&W" (32) 

and that the error behaves as 

^ k = ^ k + 0(h 2m ). 
Substituting everything back into (30), we get 

Ua = Y,Ps V s k *s+°(h 2m )- (33) 

s 

Taking together (25), (27) and (33), one gets: 

U = Uf + 0(h 2m ); I/fE^pJ^fJ (34) 

s 

where "/" stands for "filters". 

As explained in Appendix A, in the case of the Daubechies scaling functions the 
quadrature filter wi and the reconstruction filter W\ are identical. Then (34) assumes 
the following simple form: 

t^f^EW; = (35) 

This formula can be computed very fast: when c k = s k , the number of flops per 
grid point is just 4m + 3. One can use it in three dimensions too. In that case the 
quadrature and reconstruction filters are tensor products of the one-dimensional ones. 
Thus one only needs three convolutions with filters of length 2m per grid point for 
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the calculation of the potential energy, thus 12m + 3 flops. This is clearly better than 
the 32m 4 + . . . result for the triple product method. 

The formulas (34), (35) have a global nature; they characterize the approximation 
error over the whole domain of the wave function. In the next two subsections we will 
choose a local point of view instead, and prove that Eqs. (34), (35) hold for unbounded 
V(x) too. 

4-2 The matrix elements of the potential energy in the case of the 
Daubechies family. 

If we substitute the wavefunction expansion (12), (18) into the energy expression (17), 
we get 

U = J2 < J i4'*f U tj = [ dx<f>^x)V(x)^(x). (36) 

Similarly, the approximate energy (35) in the case of the Daubechies family has 
the form 

U t = £ c*a}0f ■; U\ 3 = £ Vfwt-m-j (37) 

where wi is the quadrature filter for the Daubechies family. 

It is enough to consider the matrix elements Uo q ; UL since others can be obtained 
by shifting the potential. Let us assume that V(x) = x l . Then, after going to the 
variable y = x/h under the integral one obtains: 

U 0g = J ^(x^^dx = h l K qt - K qt = J dy^^y 1 . (38) 

Similarly, one can check that 

Ul^h'K^ (39) 

i 

Now, if the potential V(x) is a smooth function, one can expand it into a Taylor series 
around the origin: 

2m- 1 y(p)(Q) 

V{x)= V ^—pL x p + 0{x 2m ). (40) 
P =o pl 

Our aim is to use wavelets for the electronic structure calculations with pseu- 
dopotentials, and the local part of the Gaussian pseudopotentials introduced in [29] 
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is smooth (has infinite number of continuous derivatives). The smoothness of the po- 
tential is best exploited if the grid is fine enough. If necessary, one can also increase 
the accuracy of our approximation without changing the number of basis functions, 
by going to the doubly dense real space grid. This is discussed in more detail in the 
very end of subsection 5.1. The multiscale nature of wavelets allows one to adapt the 
grid resolution locally, which will be described in section 5. 

Note however that one needs only the grid values of the potential for the actual 
calculations (according to Eq. (35)). The Taylor expansion of the potential (40) is 
presented here only to analyze the errors. 

Plugging the expansion (40) into (36) and taking into account (38), we get: 

2m- 1 y(p)(o) 

U °« = E -^r-h?K qp + 0(h 2m )- (41) 
p=o P- 

In the same way, plugging (40) into the right equation of (37) and using (39), one 
obtains 

ulq Jf^Y^ hPK ^ + 0{h2my 

P =o P- 

It is in the coefficients K qp that the approximation differs from the exact calcu- 
lation. Eq. (39) can be considered as a quadrature approximation of the integral 
(38). 

It would be convenient if the values of K qp and coincided. It would explain the 
smallness of error in (33), since in that case, C/^ and would differ only by 0{h 2m ). 
However, for the cases we checked (Daubechies-6 to 18) the values of K turned out 
to be different. Thus, the smoothness of the potential alone does not explain the 
smallness of error (34). In the next subsection we will see that one needs also the 
smoothness of the wavefunction to explain that. 

Note that the strategy used in [14] was to find an approximation of the matrix (36) 
such that the corresponding K v q be exactly equal to (38), i.e., that the approximation 
is exact for the polynomial potentials. This is a flexible scheme but it requires more 
computer resources than ours. The reason is that in the scheme of [14] one has to 
apply two-index filters at each point. Therefore we have the same problem as with 
the triple product formula (20): the number of flops per grid point is iV 2 in one 
dimension and iV 4 in three dimensions, where N is the dimension of the matrix 
which plays the role of (21). 
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4-3 The gradient of the potential energy in the case of the Daubechies 
family. 

The gradients of the exact and approximate energies (36) and (37) w.r.t. cf are 

where the matrix elements are discussed in the previous subsection. In the case of 
the potential energy functional, with the coefficients cf instead of in (36) and (37), 
there should be an additional factor of 2 in (42). 
From (34) one could guess that 

dU dUf 



+ 0(h 2m+1/2 ) (43) 



since there are no preferred points in space for the energy expressions (36) and (35), 
and thus the energy error should be "smeared" smoothly over the grid points. Later 
in this subsection we will see that (43) is indeed satisfied. 

The condition (43) is very important. It means that if we minimize U{ using some 
numerical method (steepest descent etc.) then the gradients along which we change 
the wavefunctions will be very close for the exact and approximated energies, and 
thus the results of minimization, the two ground state wavefunctions, will also be 
close. 

On the other hand, Eq. (34) that was derived above only for a bounded V(x), 
follows from (43). What is more, it follows from (43) that (34) is satisfied even if we 
do not require the function <3>(x) to be smooth. 

At first we will prove (43) for the partial case of i — 0, \l/(x) = x l and V(x) = x f , 
from which we then will easily make extension onto the general case. With the above 
assumptions, acting similarly to the proof of (100) one can check that 



J = h l+1 ' 2 £ C?fM l _ u (44) 



M=0 



where C" are the binomial coefficients, and M s are the scaling function moments (8). 
Substituting (38) and (44) into (42), one obtains: 

RTF 1 

™ = h l+t + l/2 AH . Au = £ C u Mi _ u £ K . tj u (45) 

° C u=0 j 

Similarly, in the case of the approximate energy, 

f)TT I 

ii = h^Ai; A\ t ^ £ CtM^ E K%f. (46) 
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The coefficients (45), (46) satisfy the following equalities: 

A lt = M l+t , l,t<m; A l0 = M h I < 2m; (47) 
A\ t = Mi +t , t + m<2m (48) 

where M s are the scaling function moments (8). Their analytical proof can be found 
in Appendix B. 

Since we were unable to find an analytical proof of (47) for the remaining relevant 
values of indices: 

Au = M l+t , m<t<2m (49) 

we checked it numerically with Mathematica for the Daubechies wavelets of the orders 
from 6 to 18. We checked both the extremal phase and least asymmetric Daubechies 
families. The relative error of (49) can be made arbitrarily small by increasing the 
precision of wavelet filters. We went down to 1CT 40 . 
Thus we see that for the Daubechies wavelets, 

A lt = 4 = M l+t ; || = g = h l+t ^ 2 M l+t ; I + 1 < 2m (50) 

for arbitrary polynomial V(x) and ty(x) with the sum of their degrees less than 2m. 
The second equality follows from (45), (46). 

In the remaining part of this subsection we will infer (43) from (50). Let us 
consider the exact gradient. It is given by the left formula (42). The matrix element 
U j has been calculated in the previous subsection. Now let us consider the scaling 
function expansion coefficient Sj (22) for a smooth function ty(x). Since the latter is 
smooth, it can be expanded into the Taylor series at the origin: 

y( x ) = 2 \2^^1 x p + 0(x 2m ). (51) 
P =o P- 

Plugging (51) into (22) and using (44), we get: 



2m- 1 \U{p)(AU\ p 

s * = £ ^_^H h P+i/2 £ C^Mt + 0{h 2m+1 ' 2 ). (52) 

p=0 P- t=0 

Substituting (52) and (41) into (42), we get: 

orr 2m-l Us+1/2 s 

™ = £ C P ^\0)V^- p \0)A PtS _ p + 0(h 2m+1 / 2 ) = 

OC s =0 31 p=0 

2m- 1 Us+l/2 

= —T~ G {S) (0)M S + O(h 2m+1/2 ) = g k G + 0(h 2m+l/2 ) (53) 
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where we have used (50). As before, G(x) = V(x)^(x), and g\ = / G(x)(j) k (x)dx. 
The analogue of (53) for an arbitrary gradient component ^ can be obtained by 

replacing all the by ih. Then our error estimate for the potential energy of Eq. (27) 
follows from the formula 

i UL -i 

For the gradient of the approximate energy one can exactly repeat the above steps, 
just using Kj t instead of Kj t (since one starts from (39) instead of (41)). The result 
will be identical to (53); thus, (43) is proved. 



5. The case of adaptive resolution. 

In the actual calculations, most of the wave function coefficients for wavelets on fine 
levels are very small. An advantage of wavelets is that the coefficients smaller than 
certain threshold can be set to zero [6], [7]; this is called dynamical adaptivity. We 
will use here the so called static adaptivity when one constrains the fine wavelets to 
be zero at the some parts of configuration space - e.g., far from the atomic cores. This 
effectively leads to grid resolution slowly varying in real space. Since our quadrature 
approximation involves some extra care of the boundaries between regions of different 
grid resolution (see below), it is more suited for the static adaptivity where one can 
make sure that these boundaries have simple form. 

The simplest two-level example will be considered below. 



5. 1 The energy error for the adaptive approximation. 

In this Section we will extensively use the bra and ket notation for the wavefunctions 
(12),(18): 

|c fc >=£cM>; (55) 

i i 

The ket vectors belong to the space Vk (2). It is finite-dimensional after we apply 
boundary conditions, whether periodic or non-periodic. 
The energy expressions (36), (37) then can be written as 

U =< c k \U\s k >; U t =< c k \U^\s k > . (56) 

Note the difference between the exact and approximate operators. U is the operator 
of multiplication by V(x) in a Hilbert space. On the other hand, U k is an operator 



17 



in the finite-dimensional space taking into account the boundary conditions. Its 
matrix elements are given by (37). 

In the present Section we will consider two resolution levels simultaneously. In 
addition to the resolution level k, considered up to now, we have a more coarse level 
k — 1 together with the vectors associated with it: 

\c k - 1 >= E^Mtf" 1 >; Id*" 1 >=E d ? _1 l^ _1 > • ( 57 ) 

i i 

For some time we will discuss only the vectors corresponding to &i(x); those re- 
lated to ^i(x) satisfy the same relations with the symbols c, d replaced by s, b where 
appropriate. 

Similarly to (22), the coefficients in (55), (57) can be expressed in terms of the 
smooth function $(x): 

c\ = J tf(x)*(x)dx- t c*" 1 = J <f>t\x)<S>(x)dx; 

dt 1 = J ^~ 1 {x)^{x)dx. 

The bra and ket notation is convenient because in accordance with (5), the forward 
and backward wavelet transformations can be written as 

\c k >= Ic*" 1 > +\d k - 1 > . (58) 

The left and right parts of the above equality differ only in the choice of basis func- 
tions: {<$} at the left and {<^~\ V^ -1 } at th e right. 

The adaptive approximation consists of replacing the second relation in (57) by 

i^^e^-v-- 1 ^ (59) 

In this paper we consider the so-called static adaptivity, where one chooses some fixed 
predetermined region D of space such that outside D the wavelet coefficients d\ are 
known to be small and can therefore be neglected. We will call D the fine region. 
Replacing > in (58) by (59), one can form the adaptive vector 

\c k >= (c^ 1 > +\d k - 1 > . (60) 

In other words, if we forward transform \c k >, then discard the wavelet part outside 
the domain D according to Eq. (59) and then backward transform the result, then 
we get \c k >. 

The discarded wavelet part can also be written as a vector: 



(61) 
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so that, taking into account (59), 

Id*- 1 >= Id*" 1 > +\d k ~ 1 > . (62) 

The exact energy expression (56) in the adaptive case assumes the form 

U & =< c k \U\s k > (63) 

where U is the matrix (36), and \s k > is formed from \s k > in the same way as \c k > 
is from \c k > (by discarding the wavelet part outside the same region D). For the 
approximate energy we make a similar ansatz: 

[/ f a =< c k \t) k \s k > (64) 

where the matrix elements of U k are given by (37). 

Let us estimate the error of the energy expression (64), i.e., its difference from (63). 
It follows from (60) and (62) that 

\c k >= \ c k > -\d k - x > (65) 

and, similarly, for \s k >: 

\S k >= \ s k > -Ifo*- 1 > . (66) 

Plugging (65) and (66) into (63), (64), we get 

= < c k \U\s k > - < d k - 1 \U\s k > - < c k \U\b k - 1 > + < d k - 1 \U\b k - 1 >; 
= < c k \U k \s k > - < d k - 1 \U k \s k > - < c k \U^\b k ^ > + 

+ < d k - 1 \u k \b k - 1 > . 

Therefore, the error of the adaptive energy approximation has the form 

SU* = U a - = 5Ux + 6U 2 + 5U 3 + 5U 4 ; 8U X =< c k \U\s k > - < c k \Uj*\s k >; 
5U 2 = < d k - l \U k \s k > - < d k ' 1 \U\s k >; 5U 3 =< c^U^- 1 > - < c k \U\b k - 1 >; 
5U± = < d k ~ 1 \U\b k - 1 > - < d k - 1 \U?\b k - 1 > . (67) 

One can show that all four error terms are small: 

5U a = 0(h 2m ). (68) 

First, 5Ui = 0(h 2m ) because this is the error term in the non-adaptive case. Then, 
let us consider the term 

su 2 =< d^Kuf - u)\s k > . 
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One can show that it contains the non-adaptive gradients: 

Fit J 

U k \s k >= U k £ s)\4>) >= £ \4> k > U{ 3 s) = £ \<j> k > 1 



dc k 

and the same for U\s k >. Therefore, 

(^-f»i^>=i:i*f>(§-||)=o(" 2m ) 

because of (43). Also, the wavelet coefficients asymptotically behave as 

d k ~ l = £>(/i m+1/2 ), so Id*" 1 >= 0{h m ) (69) 

(see [30]). Therefore, 8U2 = 0(h 3m ) and this term can be neglected. Since the 
potential energy matrices (36), (37) are Hermithean, the same applies to the third 
term. 

Eq. (69) also determines the asymptotic behavior of the fourth term: 8U4 = 
0(h 2m ). The fourth error term is formed by wavelet coefficients on the k — 1 level, 
so it is the main source of error. However, it vanishes outside D. On the other 
hand, inside D one can use the non-adaptive error estimate (34). Thus the main 
contribution to the energy error comes from the boundary of D where neither of the 
above arguments applies. 

An interesting partial case is when the fine region D is empty. Then, all the 
wavelets on the finest level k are discarded, so the wavelet vector (59) is zero. Thus 
Eq. (60) is simplified to \c k >= |c fc_1 >, which is equivalent to using just the scaling 
functions on the level k — 1. However, we still use the quadrature (10) on the level 
k, and evaluate the wave function on the grid with constant h = 2~ k , which is twice 
denser than the grid on which the scaling functions are defined. This is useful if we 
want to determine the potential energy more precisely without raising the number of 
the degrees of freedom in the variational minimization. 

5.2 The matrix elements for the adaptive case. 
If one substitutes (65), (66) directly into (63), one obtains: 

£/ a = < c k - 1 \U\s k - 1 > + < d k - 1 \U\s k - 1 > + < c k - 1 \U\b k - 1 > + 

+ < d h - 1 \TJ\~b k - 1 >; (70) 

U? = < c k ~ 1 \U^\s k - 1 > + < d k ~ 1 \U k \s k - 1 > + < c k - 1 \U f k \b k - 1 > + 

+ < d k - 1 \U k \b k - 1 > (71) 
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with the elements of the exact matrix having the form 

u% =< ^U^- 1 >= J <p k i -\x)V{x)<p k r\x)dx ] 
u cd ^< ^\(J\^ > ; u dd ^< > . 

The approximate matrix elements are a bit more difficult to derive: using the 
backward transformation and (37), we get: 

< = < tf-'WM' 1 >= E V t k v t -2iV t - 2j = J2 V£ +s V s Vs-2 g (72) 

t s 

where q — j — i, and a new filter is defined: 

v s = J2 hw s -i- (73) 

Suppose that V(x) = x l . Then, we have (38) and an analog of (39): 

u% q = (2h) t K qt ; vg = (2h)%; k{ t = ]T( S /2)W 2g . (74) 

i 

For a general potential, we have (41) 

2m- 1 y(p)(r)\ 

< q = E —ym p K qp + <D(tf m ). (75) 
p=o P- 

For other matrix elements the same formula can be used, just with wavelets instead 
of scaling functions in the definition (38) of K qp . 
If we use a variant of the Taylor series (40): 

V s k = V(sh) = ^ X ^f^{2hY{s/2f + 0{h 2m ) (76) 
p=o P- 

and combine it with (74), we will obtain an analog of (75): 

< = 2 ^^^{2h)% p + 0{h^). (77) 
p=o P- 

The matrix elements containing wavelets have the same form, the only difference 
being that the filters gi are used instead of hi where appropriate. 
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5.3 The gradient error for the adaptive approximation. 
The derivatives of (70), (71) w.r.t. cf -1 and d k ~ l have the form: 

J£r = + = - £ <r ; (^) 

-2^1 = V ug^" 1 + V w^" 1 = - V u*^" 1 (79) 

where C/ =< c k \U\s k > is the energy with all wavelets kept, i.e., for the non-adaptive 
case. As in the non-adaptive case, the expressions for the gradient of a quadratic 
form would differ just by a factor of 2. 

For the approximate energy, the gradient expression is similar. The matrix el- 
ements for the exact and approximate energy have the form (75), (77), respectively. 
Since the filters k v q have finite and short length, the gradient expressions (78), (79) 
are nearly local. Therefore, one can consider their features depending on the region 
where c k ~ l is located. 

In the fine region D (see (59) and the text after it) the b k ~ x wavelets (61) are zero. 
Therefore, Eqs. (78), (79) reduce to their non-adaptive counterparts. If we apply the 
backward wavelet transform to them, the result will be (42). Then it follows from 
(43) that the difference between the gradients of the exact and approximated energy 
will behave like 0(h 2m ). Also, all the arguments of subsection 4.3. apply. 

Now there is the region far from D where all the wavelets are discarded. We will 
also call it "coarse region". The wavelets (59): d* -1 ,^ -1 are zero there. Thus, in the 
coarse region (78) assumes the form 

= <*r- (so) 

The condition (79) is not applicable in the coarse region, since the c^ -1 wavelets are 
zero there, and the potential energy does not depend on them. 
We will prove below that in the coarse region, 

r)TJ a r)TJ a 

|L^ + 0(/ ^, (81) 

In the partial case when D = discussed in the end of subsection 5.1, the above 
equation plays the role of (43). 

In the general case there is also the border region between the fine and coarse ones. 
The error estimates (43), (81) were derived under the assumptions that either the d k 
or d k wavelets are zero, correspondingly. Since in the border region neither of these 
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assumptions is guaranteed, the estimates (43), (81) may no longer hold there (and 
the numerical tests imply that they do not hold indeed). Thus in the border region, 
one has just (78), (79), with all terms being nonzero. From the right equations of 
(43), (81) one can conclude that the error at the border is bounded by the magnitude 



of bj coefficients which behave according to (69). The other, non-adaptive-like part 
of error behaves like 0(h 2m ), because of (43). Therefore, 



dc\~ v dcf ' " v " ddt 1 ddf 

in the border region, so it is the main source of error. The small energy error (68) is 
consistent with the above because the rapidly oscillating part of the gradient at the 
boundary region has small smooth component and gets drastically diminished when 
multiplied by a smooth function, according to (54). 

In the remaining part of this subsection we will prove (81) in the coarse region. 
The proof is completely analogous to that of (43), so we will not go into much detail. 
First, let V(x) = x l and ty(x) = x l . Then, it follows from (50) that 

° U \ = (2h) t+l+1 / 2 M l+t . 



del 

Let us prove that the same is true for f/ f a . Substituting the potential energy matrix 
from (74) and (52) for the level k — 1 into (80) we obtain the following analog of (46): 



uc o t=o q 



Similarly to (50), one can show that 

f)U a f)TJ a 

a\ t = M l+t] = ^ = {2h) l+t+1 ' 2 M l+t] l + t<2m. (82) 

For the values / < m, Eq. (82) can be proved along the lines of Appendix B. Since we 
were unable to find an analytical proof of (47) for the remaining values of /, we checked 
it numerically with Mathematica for the Daubechies wavelets of the orders from 6 
to 18. We checked both the extremal phase and least asymmetric families. Same as 
with (49), the relative error of (82) can be made arbitrarily small by increasing the 
precision of wavelet filters. We went down to 10~ 40 . 

Combining (82) with the Taylor expansions (76) and (51) for the level k — 1, we 
get an analogue of (53): 



f)T[a 2m- 1 (r)l 1 \s+l/2 s 

p=0 
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2m- 1 /ol\s+1/2 

= £ GW (0)M S + 0{h 2m+1 / 2 ) = g^ 1 + 0{h 2m+1 ' 2 ) = 

= |^ + 0(^ +1/2 ). (83) 
Eq. (81) is thus proved. 

6. Reducing the computational cost in the coarse 
region. 

6. 1 The adaptive quadrature for the product of a scaling function and 
a smooth function. 

The approximation derived above is not numerically efficient in the coarse region. 
Namely, the approximate energy matrix elements are defined by (56) there. However, 
in the non- adaptive approximation for the level k — 1, they would be given by (72) 
, which requires roughly twice smaller number of calculations because the filter uj\ is 
shorter than v\ The asymptotic behavior of the energy error of these two approxima- 
tions is the same: 0{h 2m ). Therefore, it would be desirable to use (72) in the coarse 
region. 

This situation is analogous to applying the quadrature (10) on the level Ho a 
two-scale function G(x) that has "fine region" where one should use the quadrature 
(10) on the level k and "coarse region" where one should use it on the level k — 1. 

Thus we have an adaptive quadrature. Recall that the adaptive approximation 
considered in the previous Sections consists of setting the wavelet coefficients in the 
coarse region to zero. In addition to that, the adaptive quadrature is defined on a 
coarser grid in the coarse region. 

The question is - what quadrature should we use at the boundary of the coarse 
and fine regions? Our recipe is to use the fine quadrature, along with backward 
transformation, with the filter (73): 

/ ' G(x)<ftr\x)dx = VhJ2viGt 2r + 0(h M+3 / 2 ) = VhJ2v s -2rG k s + 

J I s 

+ 0(h M+3 / 2 ) 

and go to the coarse quadrature only at some distance from the boundary. 

An alternative (and maybe more rigorous) recipe [14] is to use a non-uniform 
quadrature at the boundary. However, 
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• The non-uniform quadrature is difficult to program since the points in space 
near the boundary should be dealt with individually. 

• Our error estimates for the potential energy quadrature hold only for the uni- 
form case and it is not trivial to extend them onto the case of non-uniform 
quadrature of [14]. 



6.2 The energy gradient in the coarse region. 

At present, the approximate energy gradient in the coarse region is given by 

dU a 



ccf k— 1 



Q k-1 Z-^i ij 3 

where the matrix elements are given by (72). In general , the gradient is given by 
(78),(79). 

The energy gradient on the k — 1 level in the non-adaptive scheme is given by (42) 
with the matrix elements (37) (for the level k — 1): 



dc 



k-i 



One could define a "quasigradient" vector vf 1 with the following components: 



9Ut _ v Q k ~ l TT lk ~ l i d D' 



where the set D' is D plus the points that are no more than a distance a (empirically, 
a = 3mh is enough) away from D. The word "quasigradient" means that the vector 
(85) is assembled from gradients (78), (84), but is not necessarily equal to the gradient 
of any function at all. 

We can use now use this vector instead of the gradient in the iterative minimization 
algorithms. 

The quasigradient vector (85) is a good approximation of the adaptive gradient 
(78). For the points % e D' their components coincide. For the points % ^ D', they 
coincide if V(x) and ty(x) are polynomials with the sum of their degrees smaller than 
2m, because in that case both the quasigradient and the adaptive gradient coincide 
with the exact gradient (see (50), (82)). 

For the general potentials and wave functions, the gradient (78) assumes the form 
(83). On the other hand, the quasigradient will have the form (53) for the level k — 1, 
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which is identical to (83). Therefore, the quasigradient components are asymptotically 
close to those of the adaptive gradient (78) and its approximate version: 

H - v\ = 0{h^l% % ~ = 0(/* 2m+1/2 ). (86) 

Now let us consider the energy. Eq. (64) can be rewritten in the form similar to 
(54): 



Oct 1 ' V"* ddt v 

One can define a likewise energy expression based on (85) 

^ = E^- 1 + E*" 1 ^r ( 87 ) 

i i 

where "e" stands fort for "efficient". Then it follows from (86) that 

XJl = + 0{h 2m )- U't = U a + 0{h 2m ). 

Thus, the energy expression (87) is also a good approximation of (64), (63). 

Note that the gradient of (87) does not coincide with the quasigradient (85). They 
are different only at the boundary of D'\ the former behaves badly there, while the 
latter is smooth. This means, in particular, that simply minimizing (87) would result 
in a wavefunction that is not smooth at the boundary of D' , but oscillates there 
rapidly. 

The quasigradient can be seen as a projection of the gradient of (87) onto the 
space of functions that are smooth at the boundary of D' . One can then assume that 
we minimize the energy (87) in the space of such functions. This remedies the above 
problem of singularity at the boundary. 

Alternatively, one can consider the quasigradient as an approximation of the gra- 
dient (78) and (87) as an approximation of (64). We have shown above that the 
error of such approximation is asymptotically small for smooth wavefunctions and 
potentials. 

7. Charge density and products of functions. 

In the density functional theory one needs to express the charge density 

p{x) = (tf j(ar) 
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as a linear combination of scaling functions. In the Hartree-Fock method, one uses 
products of functions of the form (12): 

^ I {x)^ I {x) = Y.c k s)<p k l {x)<P%x) (88) 

and it would also be useful to have that product expanded in terms of scaling functions 
or wavelets, not in terms of scaling function products: 

$ / (^)^/(^) = E/M'^) + ^) 
j 

where SF(x) is small. In general, one can expand the product in scaling functions of 
the order different from m (that for <&i(x), or even from some other wavelet 

family. 

The obvious way to obtain the coefficients f!f is to expand the left part in scaling 
functions according to (6). The coefficients can then be obtained by convolution with 
the matrix of triple products, similarly to (20). However, then we would have the 
same problem as with (20): too many flops per grid point, N 2 in one dimension and 
N 4 in three dimensions, where N is the dimension of the triple product matrix. 

Fortunately, we can suggest an alternative way to get the density: 



7.1. Using the average grid values of &(x) and *ff(x). 
Let us define the following coefficients: 

F k = = $(^)^) + o(h 2m ) 

where the average values are defined in (32). 

Then one can approximate the product (88) with the function 

F s (x) = h^2Fj c S(x-ih). (89) 

One can easily show that this function reproduces the multipole moments of the 
product (88): 

f F s (x)x t dx = hJ2 ^{jh)^{jh)f = f <!> I (x)x t ^ I (x)dx + 0(h 2m ) (90) 
J . J 

for all integer t > 0. The second equality is a partial case of (34). 

What is more, if ty(x) is a polynomial of the degree / and I + 1 < 2m, then (90) 
is satisfied exactly. This is a consequence of (50). 
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However, the function (89) is not smooth, in contrast to the product (88) that it 
should approximate. To remedy this, one can go from (89) to 



where 7 (y) is the interpolating scaling function of the order L. Its first L moments 
are the same as those of a delta function [9], [44]. Therefore, Fj(x) satisfies (90) too 
(fort<L). 

On the other hand, if L > m, then 



so the approximated density is close to (88) at any point. 

This approximation is much faster than that the triple product calculation. How- 
ever, the moment conservation (90) is not exact. To make at least the total charge 
of a single electron equal to one in this scheme, one has to scale the coefficients F k 
accordingly. One can also make a forward wavelet transformation to the scaling func- 
tion coefficients on some coarser level and scale only them; this is enough to get the 
correct total charge. 

8. Application to the harmonic oscillator. 

8.1. The non-adaptive case. 

In this Section we will test the above approximations by finding the ground state of 
a unit mass and frequency harmonic oscillator in one dimension. We use nonperiodic 
boundary conditions where the wavefunction was set to zero outside the interval [- 
16:16]. The variational ground state is obtained from Eq. (15), with the definitions 
(13), (11) and V(x) = x 2 /2. Since for the oscillator potential, (19) is exact, we can 
use the triple product method (20) without invoking additional errors. 
Our approximation for the potential energy has the form (35): 



= ^(x/h - i) 



F T (x) = $(x)*(x) + 0{h L ) = ^(aOttjfc) + 0{h m ) 




Accordingly one can define the approximate RR functional 



T{{4)) + U t ({c 



}). 



^ = mini? f ({ Ci fc }); 



(91) 



< c k \c k > 



> = argmini^lc^}). 
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We performed minimization of (13) in the exact and approximate cases for least 
asymmetric Daubechies - 2m with 3 < m < 8. We used the steepest descent method 
with diagonal preconditioning and gradient feedback. The results for Daubechies-8 
and Daubechies-16 are shown on Figs 2-5. On the x axis we have the inverse grid 
interval: h^ 1 = 2 k . 

In Figs 2,3, on the y axis we have the deviation of the variational ground state 
energy Eq (15) from the exact result (which is 0.5 for the unit oscillator) and the 
deviation of the approximate ground state energy Eq (91) from 0.5. 

Also shown is the difference of the variational and approximate ground state en- 
ergies (91) and (15): 

$E app = Eq — Eq 

which we call the approximation error. The graphs are on a double logarithmic scale, 
and they have a distinct linear part. The slope of the linear part is consistent with 
(16). For the lower orders, the approximation error is one or two orders of magnitude 
smaller than (16). Its slope in the linear region is equal to 2m, suggesting that 

5E app = 0(h 2m ). (92) 

This behavior is similar to (34). However, note that here we compare the energies for 
different (although close) quantum states \cq > and \c k > defined in Eqs. (15) and 
(91), while in the Section 3 we computed the exact and approximate energy for the 
same state. 

The approximation error grows with m faster than the variational error, so that for 
Daubechies-16 they become of the same magnitude, and have the same slope 2m — 2. 
In Figs 4,5, on the y axis we have the following quantities: 



sc k = ^£(4) - c k g y- 5c k app = ^J2(c k - c^y. (93) 

We will call the first quantity the variational error and the second one the approxi- 
mation error. The coefficients 

c*g = C j exp(-x 2 /2)(f) k (x)dx 

correspond to the scaling function expansion of the exact ground state of the unit 
oscillator in the space (2). The factor C is chosen such that J2i{c k g ) 2 = 1. 

We could not find the evaluation of the quantities (93) neither in [42] nor in the 
previous papers describing the application of wavelets to the Schroedinger equation 
[6], [14]. The quantity of interest in [42] was something else: the norm of the total 
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difference between the exact ground state of (11) (in this Gaussian) and 

y.;< J ;o';U): 

||7r- 1 / 4 exp(- a ; 2 /2) - £c*tf (s)|| = 0(h m ). (94) 

i 

The above quantity includes the fine wavelet part (7) of the Gaussian that also behaves 
as 0(h m ). 

For lower values of m, the slope of the graphs of (93) in the linear region is 2m — 2 
and 2m correspondingly. The approximation error in the linear region is again one 
or two orders of magnitude smaller than the variational error. This suggests the 
asymptotic behavior 

5c k = 0{h 2m - 2 ); 5c k app = 0{h 2m ). (95) 

It follows from the first estimate (95) that the wave function error (94) is dominated 
by the fine wavelet part. The second estimate of (95) means that the minima of the 
exact and approximate RR functionals are very close both to each other and to the 
projection of the exact solution onto the space of the scaling functions at resolution 
level k. This explains in part why the energy error (92) behaves similarly to (34). 

For higher values of m the error reaches the machine precision range before the 
slopes reach their asymptotic values. Still, the approximation error decays faster than 
the variational one, and the slopes differ by 2, approximately. 

For the lower values of m there is no difference between the behavior of the least 
asymmetric and extremal phase Daubechies scaling functions. However, for m > 10 
the approximation error for the extremal phase Daubechies approaches the variational 
error. The convergence of the steepest descent iterations becomes very bad. For the 
symmetric family nothing of that happens. Thus for the bigger values of m our 
method is not applicable for the extremal phase family. This is the reason why we 
prefer the least asymmetric family in general. 

In a practical electronic structure calculation the exact minimization is too costly 
numerically, but one still needs a criterion of accuracy of our approximation. We can 
use for that purpose the gradient of the exact RR functional (13) at the minimum of 
the approximate one (91). 

8.2. The adaptive case. 

We will impose adaptivity in the following way: the minimum of the RR functional 
(13) will be sought in the class of coefficients {cf _1 ,df -1 } such that the wavelet 
coefficients are zero outside the interval [—k, k]: d k ~ l = for |2 1_fe i| > k. We choose 
k — 1.5 as an illustration. The resulting energy and wavefunction are, 

E = mmR({c k -\d k - 1 }) = mmR{{d k }); \c k >= argmini?({^}) 
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where the coefficients c\ are the backward transformation of {c\ 1 ,d\ 1 }. 
Their approximate counterparts are 

E { = min J R f ({c i fe " 1 ,^" 1 }) = min R t ({c?}); |c f fc >= argmin Rf({$}) 

where in the approximate RR functional the potential energy part is treated according 
to the Sections 5-6. 

The graphs of the energy and wavefunction error for the (least asymmetric) 
Daubechies-8 and Daubechies-16 are shown on Figs 6-9. The main differences from 
the non-adaptive case are: 

• The slope of the variational wavefunction error is now approximately m + 1/2. 
The asymptotical slope of the approximation error exceeds that of the varia- 
tional error by 2, roughly. Both wavefunction errors are localized at the bound- 
ary. 

• The approximation error both for the energy and wave function for small wavelet 
orders is an order of magnitude smaller than the variational error, now for the 
bigger values of h too. 

• The behavior of the approximation error is improved for the small k, compared 
to the non-adaptive case. 

9. Conclusion. 

In the present work we propose a quadrature for the evaluation of the potential 
energy functional when the wave function is a wavelet approximation of some smooth 
function. We used the Daubechies family but the results can be extended onto others. 
With our algorithm, the potential energy can be calculated using only one-dimensional 
convolutions and filters, in contrast to the existing methods. The resulting potential 
energy differs only insignificantly from the exact value. The algorithm is extended 
onto the case of varying spatial resolution (adaptivity). As a numerical test we 
calculated the ground state energy and wave function of the harmonic oscillator in 
1 dimension for the least asymmetric and extremal phase Daubechies wavelets with 
orders from 6 to 16. We performed the minimization of the RR functional with the 
potential energy calculated by our method and compared the resulting energy and 
wavefunction with those obtained from the fully variational minimization. In the case 
of the least asymmetric Daubechies family the approximate energy and wavefunction 
are close to the variational values. However, for the extremal phase family our method 
is reliable only for wavelet orders less than 10. 
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Our method also allows a fast calculation of a charge density for a wavefunction 
expressed in the wavelet basis. As a byproduct we derived a filter for reconstruction 
of the grid values of a function from its Daubechies-2m scaling function expansion. 
This reconstruction is exact for polynomials up to order 2m and the length of the 
filter is just 2m. 

The method can readily be generalized to 3 dimensions and it is already being used 
for 3-dimensional electronic structure calculations in the framework of the BIGDFT 
project [45] which is a subject of future publication. 
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Appendix A. Reconstruction of grid values of wave 
function. 

The aim of this subsection is the reconstruction of the value of a smooth function at 
the grid points from the scaling function expansion coefficients. For this purpose we 
need to find the finite (and shortest possible) length filters Wi such that 

= ^ E W rsc k s = ^= E W i c r- q = *frr) + °(h 2m ) ( 96 ) 

where are given by (22) and h = 2~ k . 

One can start by finding filters that satisfy (96) exactly for ty(x) — x l , I < 2m 
and k — 0: 

Y. W ^sC s \A=r l (97) 

s 

where the coefficients c s [x l ] are a generalization of the moments (8): 

c s [x l ] = J x l (j) s {x)dx. 

Eq. (97) means that, though with the usual interpolation (31) using Daubechies- 
2m scaling functions one can exactly reproduce polynomials of degree not more than 
m — 1, the desired filter W\ should allow us to reproduce polynomials of degrees up to 
2m — 1. Thus for such high-order polynomials, (96) would be satisfied exactly, so the 
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average value ^ would coincide with value ^ r (rh) of the polynomial, but not with 
the value of the scaling function expansion (31) 

*i(r/0 = £#(r/0 U k l {y)^{y)dy. 



It is enough to prove (97) at the origin only: 

S^^Wscix 1 ] (98) 

s 

since then it would follow that 

Y,W r -sCs[x l ] = Y, W -^+r[x l ]=Y,W- q C q [(x + r) l ] = 
s q q 

I I 

= E^E C?r l - U c q [x u ] = Crr l ~ u 5 u = r l (99) 

q u=0 u=0 

where we changed the indices: q = s — r. To prove (98), we will use the following 
auxiliary formula: 

c s [x l ] = J x l (j) s (x)dx — I x l '<j){x — s)dx — j \y + s) 1 (p(y)dy = 

= jj2C?y l ~ u s u <P{y)dy=Y J Clts u M l „ u (100) 



u=Q u=0 



where Cf are the binomial coefficients, and M s are the scaling function moments (8). 
Substituting (100) into the right part of (98), we get: 



*£W- a c a [x l ] = Y, W ~sT, C i sUM i-n = T, C i M i-ulu(-l) u = S l ; (101) 

s s u =0 u=0 

lu = J2W s s u . (102) 



The filter W s can thus be found by solving the systems of linear equations (101), (102) 
for 7„ and then for W s . 

Now we are ready to prove (96). It is enough to do it only at the origin: 

= -^ r Y. W -sC k s = *(0) + 0(h 2m ). (103) 
Vh s 

To prove (103) one should use a version of (98) at the level k: 

4 = 1^^ (104) 
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that can be proved in the same way as (98) itself. 

Then it remains to expand ty(x) in the Taylor series at 0: 



y( x ) = + 0(h 2m ) 

P =o P- 



and plug it into (103): 



1 1 r 2m_1 

-= E w. s c^{x)] = ^= e w- s c k s E —r 1 * + otP" 

Vh s Vh S l p=0 pi 

= 4r £ W ~s 2 £ ^P-Cslx P ] + 0{h 2m ) = *(0) + 0(/. 2m ), 

the last equality following from (104). Eq. (103) is thus proved. 

Suppose we work with Daubechies-2m scaling functions. Then, Eqs. (101) are 
satisfied if 

-i u = M U) u = 0,...,2m-l (105) 
because of the equality ((B2) from [14]): 

E C s p M p _ s M s (-iy = Ml = b v - p = 0, . . . , 2m - 1 (106) 

where M p are the moments of the lazy-m scaling function. The last equality is a 
property of the interpolating scaling functions [9], [44]. Eq. (106) is satisfied both for 
the extremal phase and least asymmetric Daubechies wavelets. 

The shortest filter with the moments (105) is (10). Thus, (34) reduces to (35). 
Note however that (106) is no longer satisfied for p > 2m because the higher moments 
of an interpolating scaling function are not zero. Therefore one cannot construct a 
filter of degree p higher than 2m — 1 that would also satisfy (101) for the powers of 
x up to the p-th. Thus one cannot improve (35) so that its error scale as 0(h s ) with 
s > 2m. 

The reconstruction scheme presented above can be generalized to find the values 
of ty(x) at arbitrary points and also to find its derivatives of arbitrary order. 

Appendix B. The analytical derivation of the A\ u 
A\ t coefficients 

In this Appendix we will present the proofs of Eqs. (47), (48). 
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At first let us prove (47). Suppose, as in subsection 4.3, that V(x) = x f and 
ty(x) = x l . Then, combining the left equalities of (42) and (38), we get, for I < m, 



f| = T,SjUoj = E s ? / <p k (x)x t ^(x)dx = J <p\x)x t+l dx = h t+l+l ' 2 M m (107) 



where we used the fact that J2j Sj[x l ](f>j(x) = x l for / < m. From (107) and (45), we 
thus get the first equality of (47). On the other hand, if t — then (38) reduces to 
Uoj = 5j and thus (42) has the form 

|| = s k [x l ] = h'+V2 Mj (108) 

where the last equality is a partial case of (44). Combining (108) and (45) we get the 
second equality of (47) 

Now let us turn to Eq. (48). Plugging (39): 

U 0q = & ^WrWr-g 
r 

into the second equality of (42), we get: 

7T¥ = E S i U 0j = h * E T ' tw r E W r-q S q W] = ^ E ^W^ihr) 1 = 

ocq ■ r q r 

= h t+l+l ' 2 E w r r t+l = h t+l+l ' 2 M t+l (109) 

r 

where we have used (99) and (105). Now, Eq. (48) follows from (46) and (109). 
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Figure 1: The least asymmetric Daubechies-8 scaling function and the cor- 
responding quadrature filter. As seen from the graph, the filter values are 
close to the scaling function values at integer points. 
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Figure 2: The energy errors for the least asymmetric Daubechies-8 scaling 
functions. The slope of the variational error is 6 and that of the approxima- 
tion error is 8. 
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Figure 3: The energy errors for the least asymmetric Daubechies-16 scaling 
functions. The slope of variational error reaches 14 but that of the approxi- 
mation error reaches only 14.7 
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Figure 4: The wave function errors for the least asymmetric Daubechies-8 
scaling functions. The slope of the variational error is 6 and that of the 
approximation error is 8. 
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Figure 5: The wave function errors for the least asymmetric Daubechies-16 
scaling functions. The slope of the variational error reaches only 11.5 and 
that of the approximation error reaches 14. 



5 




1e-12 1 1 1 — 1 1 — ■ 1 ■ 1 

1 10 

1/h 

Figure 6: The energy errors for the least asymmetric Daubechies-8 scaling 
functions with adaptivity. The slope of the variational error is 6 and that 
of the approximation error is 8. The approximation error is always smaller 
than the variatonal one. 
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Figure 7: The energy errors for the least asymmetric Daubechies-16 scaling 
functions with adaptivity. The slope of the variational error reaches 13.3 
and that of the approximation error reaches 16. The approximation error is 
always smaller than the variatonal one. 
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Figure 8: The wave function errors for the least asymmetric Daubechies-8 
scaling functions with adaptivity. The slope of the variational error is 4.5 
and that of the approximation error is 7. The approximation error is always 
smaller than the variatonal one. 
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Figure 9: The wave function errors for the least asymmetric Daubechies-16 
scaling functions with adaptivity. The slope of the variational error reaches 
9 and that of the approximation error reaches 11. The approximation error 
is always smaller than the variatonal one. 
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